import numpy as np
import random
import copy

"""
本文
https://blog.csdn.net/u010866505/article/details/77877345

该文章参考：
http://blog.csdn.net/zouxy09/article/details/8537620/
http://www.cnblogs.com/jerrylead/archive/2011/04/06/2006936.html

一组序列中的数字由两个分布生成的数字混合而成。
现在估计一下这个混合分布的统计量，如均值
"""

SIGMA = 6
EPS = 0.0001


#生成方差相同，均值不同的样本
def generate_data():
    miu1 = 20
    miu2 = 40
    N = 1000
    X = np.mat(np.zeros((N,1)))
    for i in range(N):
        temp = np.random.uniform(0,1)
        if (temp > 0.5):
            X[i] = temp * SIGMA + miu1
        else:
            X[i] = temp * SIGMA + miu2
    return X


def my_em(X):
    k = 2
    N = len(X)
    miu = np.random.rand(k,1)
    Posterior = np.mat(np.zeros((N,2)))
    dominator = 0
    numerator = 0
    for iter in range(1000):
        for i in range(N):
            dominator = 0
            for j in range(k):
                dominator = dominator + np.exp(-0.1 / (2.0 * SIGMA**2) * (X[i] - miu[j])**2)
            for j in range(k):
                numerator = np.exp(-0.1 / (2.0 * SIGMA**2) * (X[i] - miu[j])**2)
                Posterior[i, j] = numerator / dominator

        old_miu = copy.deepcopy(miu)
        for j in range(k):
            numerator = 0
            dominator = 0
            for i in range(N):
                numerator = numerator + Posterior[i, j] * X[i]
                dominator = dominator + Posterior[i, j]
            miu[j] = numerator / dominator
        # print(np.abs(miu - old_miu))
        if np.abs(miu - old_miu).sum() < EPS:
            print(f"miu:{miu}, iter:{iter}")
            break


if __name__ == "__main__":
    X = generate_data()
    # print("X:", X)
    print(np.average(X))
    my_em(X)